##################################################################################################
# "Natural Disasters, 'Partisan Retrospection,' and U.S. Presidential Elections" #################
# Boris Heersink, Jeffrey Jenkins, Michael Olson, and Brenton Peterson ###########################
# Panel Analysis #################################################################################
##################################################################################################

# after going to stata, estimating the model, generating the predicted values,
# save them, and send them back here.
  
  copart <- haven::read_dta("newdata_copart_fits.dta")
  copart$yhat <- copart$yhat-copart$yhat[copart$propertypercapl_6mo==0]
  copart$Type <- "Co-Partisan"
  
  contra <- haven::read_dta("newdata_contra_fits.dta")
  contra$yhat <- contra$yhat-contra$yhat[contra$propertypercapl_6mo==0]
  contra$Type <- "Contra-Partisan"
  
  plot_dat <- rbind(copart[,c("propertypercapl_6mo","se","yhat","Type")],
                    contra[,c("propertypercapl_6mo","se","yhat","Type")])

# what is the median values?
  
  gr_corrected <- read.csv("gr_corrected.csv")

  median <- median(gr_corrected[gr_corrected$propertypercapl.6mo>0 & !
                                  is.na(gr_corrected$propertypercapl.6mo),"propertypercapl.6mo"],na.rm=T)
  
  print(median)
  
  print(exp(median))

# do some calculations

  sum(coef(gr_aug_int_new)[c("propertypercapl.6mo","propertypercapl.6mo:copart")])*median
  
  sum(coef(gr_aug_int_new)[c("propertypercapl.6mo")])*median

# make the plot

  cairo_ps("./results/fig4.eps",width=8,height=5,fallback_resolution = 1200)
  print(
    ggplot(as.data.frame(plot_dat),aes(x=propertypercapl_6mo,y=yhat,ymax=yhat+2*se,ymin=yhat-2*se,group=Type,colour=Type))+
      geom_hline(yintercept=0,colour="black",size=1.5,linetype=2)+
      geom_point(size=3)+
      geom_linerange(size=1.5)+
      geom_vline(xintercept=median,linetype=3,size=1.5)+
      scale_colour_grey(start=0.4,end=0.7)+
      xlab("Disaster Damage")+
      ylab("Estimated Treatment Effect")+
      scale_x_continuous(breaks=c(0,log(50),log(1000),log(100000),log(10000000)),
                         labels = c("$0","$50","$1,000","$100,000","$10,000,000"))+
      theme_minimal()
  )
  dev.off()
  
  png("./results/fig4.png",width=8,height=5,units="in",res = 1200)
  print(
    ggplot(as.data.frame(plot_dat),aes(x=propertypercapl_6mo,y=yhat,ymax=yhat+2*se,ymin=yhat-2*se,group=Type,colour=Type))+
      geom_hline(yintercept=0,colour="black",size=1.5,linetype=2)+
      geom_point(size=3)+
      geom_linerange(size=1.5)+
      geom_vline(xintercept=median,linetype=3,size=1.5)+
      scale_colour_grey(start=0.4,end=0.7)+
      xlab("Disaster Damage")+
      ylab("Estimated Treatment Effect")+
      scale_x_continuous(breaks=c(0,log(50),log(1000),log(100000),log(10000000)),
                         labels = c("$0","$50","$1,000","$100,000","$10,000,000"))+
      theme_minimal()
  )
  dev.off()
  
